

close all;
clear all;

tic
tot_start=tic;
global k1a 
global k1b 
global k2a  
global k2b 
global k3 
global k4a 
global k4b 
global k5a 
global k5b 
global k6a 
global k6b 
global k7 
global eta
global dx
global dx2
global ksyn %Masha dev oct2012
global k8a %rate of pheromone binding to receptor complex                                   %Masha dev oct2012
global k8b %rate of pheromone unbinding to receptor complex                                 %Masha dev oct2012
global k9 %rate of nucleotide exchange in Cdc42D->T assisted by liganded receptor RecA %Masha dev oct2012
 
Cdc42_store   = fopen('data_Cdc42T_time_course', 'w'); 
neg_store     = fopen('data_negative',          'w'); % going to hold the time points and protein when it goes negative
mass_store    = fopen('data_mass',              'w'); % going to hold total mass
Rec_store =fopen('data_receptor_time_course','w'); %storing time evolution of the receptor complex concentration on the membrane %Masha dev oct2012
RecA_store     =fopen('data_reca_time_course',  'w'); %storing time evolution of the receptor complex bound to pheromone on the membrane %Masha dev oct2012
rec_count_time=fopen('data_reca_time_course.txt',  'w'); %Debug for reception recycling and flux calculation

%% variables for time-stepping
dt    = 0.05;   % diffusion time step in seconds
%sim   = 90*100;  % number of 90 second intervals to simulate
sim   = 5*60*60.0; %Masha debug - simulation time in (s)
display(['Simulation time is ' num2str(sim/60) '(min), with a time step ' num2str(dt) ' (s)']); %Masha debug
simdt = sim/dt; % number of dt time steps
tsave = 20 %Data is going to be saved every tsave seconds

%% variables for spatial discretization
N = 100;                            % number of spatial grid points
Cellsize = 5;                       % sphere diameter (microns)
Internalcompartment = 0.7;          % relative size (area) of internal compartment vs plasma membrane 
dx  = Cellsize*sqrt(pi)/N;          % spatial meshsize (1D)
dx2 = sqrt(Internalcompartment)*dx; % spatial meshsize for internal compartment (1D)


%% initalising diffusion mechanism ________________________________________

nRsteps = 100;      % take nRsteps reaction steps per diffusion step
dt2 = dt/nRsteps;   % reaction time-step

totcnt = 1;

%% initial concentrations
version='0';
%membrane bound species
load(['Cdc42T_end' version '.mat']);
load(['BemGEF42_end' version '.mat']);
load(['BemGEF_end' version '.mat']);
load(['BemGEFc_end' version '.mat']);
load(['Cdc42D_end' version '.mat']);
load(['GDI42c_end' version '.mat']);
load(['GDI42_end' version '.mat']);
load(['GDIc_end' version '.mat']);
load(['vSNARE_end' version '.mat']);
%load(['Rec_end' version '.mat']);  %Receptor complex on the membrane %Masha dev oct2012
%load(['RecA_end' version '.mat']);      %Receptor bound to A pheromone on the membrane %Masha dev oct2012
Rec = 3.1778*ones(N); %Receptor initial concent on the memb. [microM], given cell rad 5micron, mem thickness 10 nm, tot number 6000 per cell %Masha dev sep2012
RecA = zeros(N); %Pheromone bound receptor, init conc [microM] on the membrane.


%cytoplasmic species
load(['vSNAREic_end' version '.mat']);
load(['GDI42ic_end' version '.mat']);
load(['Cdc42Dic_end' version '.mat']);
%load(['Recic_end' version '.mat']); %Newly synthesized receptors in the internal compartment %Masha dev oct2012
Recic=0;

% well-mixed cytoplasmic species
BemGEFc = mean(mean(BemGEFc));
GDI42c  = mean(mean(GDI42c));
GDIc    = mean(mean(GDIc));

%Pheromone constant concentration profile
% ee=0:1/N:(1-1/N); 
% [X, Y]=meshgrid(ee);
% Pher=0.01*(1-sqrt((X-0.5).*(X-0.5)+(Y-0.5).*(Y-0.5))); %cone, centered in the middle of the grid.
Pher=0.01*ones(N); %uniform; 

display(['Average pheromone concentration is ' num2str(mean(mean(Pher))) ' (microMolar)']);
save Pher_end0.mat Pher;

%shift peak to a particular spot, if needed for pheromone sensing
shiftpeaketc %if shift is not needed, comment this line.

%% reaction constants

mult  = 16;

k1a = 10;       %s-1,       BemGEFc -> BemGEF
k1b = 10;       %s-1,       BemGEF -> BemGEFc

% GEF
k2a = 0.16;     %uM-1.s-1,  BemGEF + Cdc42D -> Cdc42T
k3  = 0.35;     %uM-1.s-1,  BemGEF42 + Cdc42D -> Cdc42T
%GAP
k2b = 0.63;    %s-1,       Cdc42T -> Cdc42D 

k4a = 10;       %uM-1.s-1,  BemGEF + Cdc42T -> BemGEF42
k4b = 10;       %s-1,       BemGEF42 -> BemGEF + Cdc42T
k7  = 10;       %uM-1.s-1,  BemGEFc + Cdc42T -> BemGEF42

% Cdc42D on
k5a = mult*9;      %s-1,       RDI42c -> RDI42
k6b = mult*5;      %s-1,       RDI42 -> Cdc42D + RDIc
% Cdc42D off
k6a = mult*15;      %uM-1.s-1,  Cdc42D + RDIc -> RDI42
k5b = mult*1.3;     %s-1,       RDI42 -> RDI42c

%Receptor (Complex) birth and sink rates [microM/s] %Masha dev oct2012
%ksyn = 0.0021143; % [uM/s] given 0.006*Vcell volume of yeast golgi and synthesis rate 4 molecules/s per cell. Cell radius is 5 micron.
%ksyn = 0.05586*3.1778*0.01 %given the rate of exp receptor decay without RecA and nucl.exch. reactions. %Masha nov2012 
ksyn = 0.0; %no additional synthesis

%Rec + alpha <-> RecA %Masha dev oct2012
%k8a=1.0; % binding of pheromone. Simon PNAS, 2003,100(19),10764-10769
% k8a=1.0*Pher; %pheromone dependent 2D constant
% k8b=0.01; %unbinding of pheromone. Simon PNAS, 2003,100(19),10764-10769
KD8=10*10^(-3)% KD of pheromone binding reaction Raths et al. 1988 PMID 2846561 gives 20 nM, Jenness, D.D., Burkholder, A.C., and Hartwell, L.H. (1986) Mol. Cell. Biol. 6, 318-320 give 6 nM. Take an average btwn the two, 10 nM. 
k8b=0.001 %Raths et al. 1988 PMID 2846561
k8a=(k8b/KD8)*Pher; %Keep KD constant in case individual rates need to be varied

%Cdc42D + RecA ->(k9)-> Cdc42T + RecA
%k9=1*k2a;
k9=0;
display(['Rate of Cdc42D->Cdc42T, facilitated by RecA, k9 =' num2str(k9)]);

eta0=0.01;  %eta: Vm/Vc, membrane/cytoplasm volume correction
eta =eta0;

% notes in '08_endocytic_vesicle_30Dec2010.doc'
mem_depth  = (Cellsize/2)*(1-(1/((eta0+1)^(1/3))));
cyt_mult   = ((Cellsize/2)^3)/(eta0+1);
Rnew_mult  = N/2/(pi^0.5);

% compute diffusion-coeff times dt/dx^2
Dconst = 0.0045;                % diffusion coefficient
Diff1  = Dconst * dt/dx^2;

% diffusion matrix, includes dt/dx^2 factor
e        = ones(N,1); 
Dxx      = spdiags([e, -2*e, e], -1:1, N, N);
Dxx(1,N) = 1;                   % periodic boundary conditions
Dxx(N,1) = 1;
Dxx      = Dxx*Diff1;

Mxx = kron(speye(N),Dxx);
Myy = kron(Dxx,speye(N));
Hop = speye(N^2) - Mxx - Myy;   % heat operator, I-Dxx

%% initialize vesicle mechanism ___________________________________________

% Parameters related to vesicle fusion/fission

TOP = 200;              % number of bins in the window of highest Cdc42/GTP
                        % concentration
                        
cable       = [];       % hold upto MAXcable cable co-ordinates
MAXcable    = 5;        % maximum number of cables
detatchPROB = dt/(60*2);    % (time step)/(residence time),	if time step < residence time    
attachPROB  = dt;       % time step, if time step < 1, (prob 1 if time step = 1sec)

nVInlen  = N/50;        % no. of bins (1D) used to represent an exocytic 
                        % vesicle
nVOutlen = N/100;       % no. of bins (1D) used to represent an endocytic 
                        % vesicle
periodOut = 0.6/2;        % mean period for endocytosis, sec (= 15 sec actin 
                        % patch/25 patches in cell)
periodIn = 4*periodOut; % mean period for exocytosis, sec
patch_fill = 10;        % fill level for endocytosis, used for sink_case =2
patch_fill_give_up = 24;% even if the patch doesn't fill up with cargo, 
record_vSNARE = [];     % record [v-SNARE] at time of endocytosis
record_endot = [];      % record time of endocytosis
record_lifetime = [];   % record lifetime of each sink
%Masha dev oct2012
record_Rec = [];     % record [Rec] at time of endocytosis
record_RecA = [];     % record pheromone bound receptor [RecA] at time of endocytosis

endo_red_time = 15;     % 15 s before endocytosis when spots are "red"
exo_cnt = 0;            % count exocytic events
xg = (repmat(1:N,N,1)); yg = xg';  % xg, yg used to to identify window
           
% random seed
%RandStream.setGlobalStream(RandStream('mt19937ar','seed',00)); %Produces identical trajcetories! %Masha dev nov2012
RandStream.setGlobalStream(RandStream('mt19937ar','Seed','shuffle')); %Fix: clock&date as a seed -> produces uncorrelated trajectories. %Masha dev nov2012

% Finding the initial TOP bins
Cdc42T_totm = Cdc42T+BemGEF42;                % total Cdc45/GTP membrane
                                              % bistribution

[val,Tbins] = sort(Cdc42T_totm(:),'descend'); % 'val' = sorted concentrations
                                              % 'Tbins' = sorted index    
circ_mem = Tbins(1:TOP)';

circ_ind    = [];               % will hold positions of TOP bins
ves_prob    = [];               % value 1 in TOP bins, 0 otherwise
ves_outprob = [];               % value 0 in TOP bins, 1 otherwise

Ratio = 40;         % probability of a bin endocytosing within the window
                    % are 40x higher than the probability of a endocytosing
                    % outside

% probability that endocytosis occurs inside TOP bins
%endo_m = 40/89;
endo_m = Ratio * TOP/N^2 / ((1-TOP/N^2) + Ratio*TOP/N^2)

% storage for wrapping, cells chosen to endocytose, counters, and tracking Cdc42 distribution
Fwarp = zeros(N+1); % storage for use in interpolation that needs wrap around func values
endo_cell        = struct('endo_time',[],'pick_time',[],'vSNARE',[],'Rec',[],'RecA',[],'Cdc42T',[], ...
    'Cdc42D',[], 'BemGEF42',[], 'GDI42',[], 'BemGEF',[], 'posx',[],'posy',[],'dx',[],'num',0); 

%% write initial data into files
      for i = 1:N
%         fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+Cdc42D(i,:)+GDI42(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'\n');
        fprintf(RecA_store,'%5.4f ',RecA(i,:));
        fprintf(RecA_store,'\n');
        fprintf(Rec_store,'%5.4f ',Rec(i,:));
        fprintf(Rec_store,'\n');
      end
      
      % Calclate the global mass __________________________________________
     Cdc42s = GDI42c ...                                    % Cdc42D cytoplasmic content 
     + eta*(mean(mean(GDI42+Cdc42D+Cdc42T+BemGEF42))) ...   % Cdc42 membrane content 
     + eta*(Cdc42Dic*(dx2^2)/(dx^2)) ...                    % Cdc42D on inner membrane 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane           
     + eta*sum(endo_cell.Cdc42T(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42T in vesicle                            
     + eta*sum(endo_cell.Cdc42D(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42D in vesicle 
     + eta*sum(endo_cell.GDI42(1:endo_cell.num)    .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % GDI42 in vesicle
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 in vesicle 
     
     Bem1s = BemGEFc ...                                    % cytoplasmic content scaled to outer membrane
     + eta*(mean(mean(BemGEF+BemGEF42))) ...                % Bem1 membrane content  
     + eta*sum(endo_cell.BemGEF(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Bem1 vesicle content 
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 vesicle content                                                        

     GDIs = GDIc+GDI42c ...                                 % GDI cytoplasmic content 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane 
     + eta*(mean(mean(GDI42))) ...                          % GDI on outer membrane
     + eta*sum(endo_cell.GDI42(1:endo_cell.num).* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);            % GDI42 in vesicle

     vSNAREs        = eta*mean(mean(vSNARE)) ...              % v-SNARE membrane content
     + eta*(vSNAREic*(dx2^2)/(dx^2)) ...                   % v-SNARE on inner membrane content 
     + eta*sum(endo_cell.vSNARE(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % v-SNARE vesicle content 

     Recs        = eta*mean(mean(Rec)) ...              % Receptor complex membrane content %Masha dev oct2012
     + eta*(Recic*(dx2^2)/(dx^2)) ...                   % Receptor complex on inner membrane content 
     + eta*sum(endo_cell.Rec(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 
     
     RecAs        = eta*mean(mean(RecA)) ...              % Receptor complex membrane content %Masha dev oct2012 
     + eta*sum(endo_cell.RecA(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 

      %____________________________________________________________________                                                  
     
      % writing global mass into file _____________________________________
      fprintf(mass_store,'%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n',0,Cdc42s, Bem1s, GDIs,vSNAREs,Recs,RecAs);

      % calculate the internal compartment mass ___________________________
      % ic_conc = [Cdc42Tic GDI42ic vSNARE]
      ic_conc(totcnt,1) = (Cdc42Dic);   % Cdc42D  on inner membrane 
      ic_conc(totcnt,2) = (GDI42ic);
      ic_conc(totcnt,3) = (vSNAREic);   % v-SNARE on inner membrane                                      
      %____________________________________________________________________
            
      totcnt = totcnt+1;
%%

display(['Initiation stage took ' num2str(toc(tot_start)) ' (s) is done. Starting iterations...']); %Masha debug
tcycles=[];treact=[];tcargos=[];tdiffuse=[];
exotimes=[];endotimes=[];
%Masha debug jan2013
% figure(2);
% subplot(2, 2, 1);
% clims=[0 max(max(Cdc42T+BemGEF42))];
% imagesc(Cdc42T+BemGEF42, clims);
% colorbar;
% title('Cdc42-GTP concentration, [microMolar]');
% subplot(2, 2, 2);
% clims=[0 3.3];
% imagesc(RecA,clims);
% colorbar;
% title('RecA concentration, [microMolar]');
% subplot(2, 2, 3);
% clims=[0 3.3];
% imagesc(Rec,clims);
% colorbar;
% title('Rec concentration, [microMolar]');
% subplot(2, 2, 4);   
% imagesc(Pher,clims);
% colorbar;
% title('Pheromone concentration, [microMolar]');
% tstr=['Time is ' num2str(0) ' (s) out of ' num2str(sim/(60*60)) ' (hr)'];
% ah=gca;
% axes('position',[0,0,1,1],'visible','off');
% tx = text(0.5,0.98,tstr);
% set(tx,'fontweight','bold');
% axes(ah);

for loop=1:simdt                % looping over dt time steps
    tcycle=tic;
    react_step_Euler_VSS_GDI42ic   
    treact=[treact toc(tcycle)];
    curr_t = loop*dt;
    
    % Finding the TOP bins after each react_step
    Cdc42T_totm = Cdc42T+BemGEF42;                % total Cdc45/GTP membrane
                                                  % bistribution
                                            
    [val,Tbins] = sort(Cdc42T_totm(:),'descend'); % 'val' = sorted concentrations
                                                  % 'Tbins' = sorted index    
    circ_mem = Tbins(1:TOP)';
                
    circ_ind = circ_mem;                          % index of TOP concentrations

    ves_prob = zeros(N); ves_prob(circ_ind)=1;    % 1 in TOP bins, 0 otherwise
    ves_outprob = ones(N);                        % 0 in TOP bins, 1 otherwise
    ves_outprob(circ_ind)=0;                
%     pcolor(ves_prob);pause(0.1)
    tcargo=tic;
    cargo_step_WT_GDI42ic     
    tcargos=[tcargos toc(tcargo)];
    if(mod(curr_t,0.5)==0) % progress report
        display([' Current time step: ',num2str(curr_t/60) ' (min)']);
        fig2=figure(2);
%         image(Cdc42T+Cdc42D+GDI42+BemGEF42);
%         imagesc(Cdc42T+BemGEF42, [0 max(max(Cdc42T+BemGEF42))]);
%         clims=[0 max(max(RecA))];
%
%         clf(fig2);
%         subplot(2, 2, 1);
%         clims=[0 max(max(Cdc42T+BemGEF42))];
%         imagesc(Cdc42T+BemGEF42, clims);
%         colorbar;
%         title('Cdc42-GTP concentration, [microMolar]');
%         subplot(2, 2, 2);
%         clims=[0 3.3];
%         imagesc(RecA,clims);
%         colorbar;
%         title('RecA concentration, [microMolar]');
%         subplot(2, 2, 3);
%         clims=[0 3.3];
%         imagesc(Rec,clims);
%         colorbar;
%         title('Rec concentration, [microMolar]');
%         subplot(2, 2, 4);   
%         imagesc(Pher,clims);
%         colorbar;
%         title('Pheromone concentration, [microMolar]');
%         tstr=['Time is ' num2str(curr_t) ' (s) out of ' num2str(sim/(60*60)) ' (hr)'];
%         tx = text(0.4,0.95,tstr);
%         set(tx,'fontweight','bold');
    end

    if(mod(curr_t,tsave)==0) % take a snap shot every 90 seconds
      save output.mat
      for i = 1:N
%        fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+Cdc42D(i,:)+GDI42(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'%5.4f ',Cdc42T(i,:)+BemGEF42(i,:));
        fprintf(Cdc42_store,'\n');
        fprintf(RecA_store,'%5.4f ',RecA(i,:));
        fprintf(RecA_store,'\n');
        fprintf(Rec_store,'%5.4f ',Rec(i,:));
        fprintf(Rec_store,'\n');
      end
      
      % Calclate the global mass __________________________________________
     Cdc42s = GDI42c               ...                      % Cdc42D cytoplasmic content 
     + eta*(mean(mean(GDI42+Cdc42D+Cdc42T+BemGEF42))) ...   % Cdc42 membrane content 
     + eta*(Cdc42Dic*(dx2^2)/(dx^2)) ...                    % Cdc42D on inner membrane 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane           
     + eta*sum(endo_cell.Cdc42T(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42T in vesicle                            
     + eta*sum(endo_cell.Cdc42D(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Cdc42D in vesicle 
     + eta*sum(endo_cell.GDI42(1:endo_cell.num)    .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % GDI42 in vesicle
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 in vesicle 
     
     Bem1s = BemGEFc ...                                    % cytoplasmic content scaled to outer membrane
     + eta*(mean(mean(BemGEF+BemGEF42))) ...                % Bem1 membrane content  
     + eta*sum(endo_cell.BemGEF(1:endo_cell.num)   .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2) ...     % Bem1 vesicle content 
     + eta*sum(endo_cell.BemGEF42(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);        % BemGEF42 vesicle content                                                        

     GDIs = GDIc+GDI42c ...                                 % GDI cytoplasmic content 
     + eta*( GDI42ic*(dx2^2)/(dx^2)) ...                    % GDI42 on inner membrane 
     + eta*(mean(mean(GDI42))) ...                          % GDI on outer membrane
     + eta*sum(endo_cell.GDI42(1:endo_cell.num).* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);            % GDI42 in vesicle

     vSNAREs        = eta*mean(mean(vSNARE)) ...              % v-SNARE membrane content
     + eta*(vSNAREic*(dx2^2)/(dx^2)) ...                   % v-SNARE on inner membrane content 
     + eta*sum(endo_cell.vSNARE(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % v-SNARE vesicle content 
     
     Recs        = eta*mean(mean(Rec)) ...              % Receptor complex membrane content %Masha dev sept2012
     + eta*(Recic*(dx2^2)/(dx^2)) ...                   % Receptor complex on inner membrane content 
     + eta*sum(endo_cell.Rec(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content 
    
     RecAs        = eta*mean(mean(RecA)) ...              % Receptor complex membrane content %Masha dev sept2012
     + eta*sum(endo_cell.RecA(1:endo_cell.num) .* (endo_cell.dx(1:endo_cell.num).^2))/(dx^2)/(N^2);          % Receptor complex vesicle content

      %____________________________________________________________________                                                  
     
      % writing global mass into file _____________________________________
      fprintf(mass_store,'%5.4f %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\n',curr_t,Cdc42s, Bem1s, GDIs,vSNAREs,Recs,RecAs);

      % calculate the internal compartment mass ___________________________
      % ic_conc = [Cdc42Tic GDI42ic vSNARE]
      ic_conc(totcnt,1) = (Cdc42Dic);   % Cdc42D  on inner membrane 
      ic_conc(totcnt,2) = (GDI42ic);
      ic_conc(totcnt,3) = (vSNAREic);   % v-SNARE on inner membrane                                      
            
      %____________________________________________________________________
            
      totcnt = totcnt+1;
    end
    tcycles=[tcycles toc(tcycle)];      
  end  

display(['Average time per cycle ' num2str(mean(tcycles)) ' +- ' num2str(std(tcycles)) ' (s)']); %Masha debug
display(['Average time for react. part in 1 cycle ' num2str(mean(treact)) ' +- ' num2str(std(treact)) ' (s)']); %Masha debug
display(['Average time for cargo  part in 1 cycle ' num2str(mean(tcargos)) ' +- ' num2str(std(tcargos)) ' (s)']); %Masha debug
display(['Average time for diffusion in 1 cycle ' num2str(mean(tdiffuse)) ' +- ' num2str(std(tdiffuse)) ' (s)']); %Masha debug
display(['Exocytic events happened :' num2str(size(exotimes)) ' times']); %Masha debug
display(['Endocytic events happened :' num2str(size(endotimes)) ' times']); %Masha debug
%% ________________________________________________________________________

fclose(Cdc42_store);
fclose(neg_store);
fclose(mass_store);
fclose(Rec_store); %Masha dev oct2012

newv=num2str(str2num(version)+1); %after each restart, the end concentration profile will be saved into a separate file.
save(['vSNAREic_end' newv '.mat'], 'vSNAREic');
save(['GDI42ic_end' newv '.mat'], 'GDI42ic');
save(['Cdc42Dic_end' newv '.mat'], 'Cdc42Dic');
save(['Recic_end' newv '.mat'], 'Recic'); %Save last snapshot of Receptor complex in the internal compartment %Masha dev oct2012

save(['vSNARE_end' newv '.mat'], 'vSNARE');
save(['Rec_end' newv '.mat'], 'Rec'); %Save last snapshot of Receptor complex on outer membrane %Masha dev oct2012
save(['RecA_end' newv '.mat'], 'RecA'); %Save last snapshot of Receptor bound to pheromone on outer membrane %Masha dev oct2012
save(['GDI42_end' newv '.mat'], 'GDI42');
save(['Cdc42D_end' newv '.mat'], 'Cdc42D');
save(['Cdc42T_end' newv '.mat'], 'Cdc42T');
save(['BemGEF_end' newv '.mat'], 'BemGEF');
save(['BemGEF42_end' newv '.mat'], 'BemGEF42');

save(['GDIc_end' newv '.mat'], 'GDIc');
save(['GDI42c_end' newv '.mat'], 'GDI42c');
save(['BemGEFc_end' newv '.mat'], 'BemGEFc');
display(['Total time = ' num2str(toc(tot_start)/60) ' (min).']);
%time_whole = toc(tot_start)
save output.mat
